CCA_3_1_1_batch_correction_seurat_5000_23_final_4_for_html.R

## loading libraries
library(ggplot2)
library(ggpubr)

## Loading required package: magrittr

library(Seurat)
library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggpubr':
## 
##     get_legend

library(readxl)
library(corto)
library(Matrix)
library(EnhancedVolcano)

## Loading required package: ggrepel

## Initilization 
set.seed(0)
label_size=8
point_size=1
folder="~/corona_project/plots_intigrate/"
cell_type_genes=list("Bowman’s glands"=c("SOX9", "SOX10", "GPX3"),
                     "olfactory HBCs"=c("TP63", "KRT5", "CXCL14", "SOX2", "MEG3"),
                     "olfactory ensheathing glia"=c("S100B", "PLP1", "PMP2","MPZ", "ALX3"),
                     "olfactory microvillar cells"=c("ASCL3", "CFTR", "HEPACAM2", "FOXL1"),
                     "immature neurons"=c("GNG8", "OLIG2", "EBF2", "LHX2", "CBX8"),
                     "mature neurons"=c("GNG13", "EBF2", "CBX8", "RTP1"),
                     "GBCs"=c("HES6", "ASCL1", "CXCR4", "SOX2", "EZH2", "NEUROD1", "NEUROG1"),
                     "sustentacular cells"=c("CYP2A13", "CYP2J2", "GPX6", "ERMN", "SOX2"))
cell_type_genes_temp=list("Bowman’s glands"=c("SOX9", "GPX3"),
                          "olfactory HBCs"=c("CXCL14", "MEG3"),
                          "olfactory ensheathing glia"=c("S100B", "PLP1"),
                          "olfactory microvillar cells"=c( "ASCL3", "HEPACAM2"),
                          "immature neurons"=c("GNG8","OLIG2"),
                          "mature neurons"=c("GNG13","RTP1"),
                          "GBCs"=c("CXCR4", "SOX2"),
                          "sustentacular cells"=c("CYP2A13", "ERMN"))
marker_genes=c("SOX9", "SOX10", "GPX3",
               "TP63", "KRT5", "CXCL14", "SOX2", "MEG3",
               "S100B", "PLP1", "PMP2","MPZ", "ALX3",
               "ASCL3", "CFTR", "HEPACAM2", "FOXL1",
               "GNG8", "OLIG2", "EBF2", "LHX2", "CBX8",
               "GNG13", "EBF2", "CBX8", "RTP1",
               "HES6", "ASCL1", "CXCR4", "SOX2", "EZH2", "NEUROD1", "NEUROG1",
               "CYP2A13", "CYP2J2", "GPX6", "ERMN", "SOX2")
marker_genes_temp=c("SOX9", "GPX3","CXCL14", "MEG3","S100B", "PLP1", "ASCL3", "HEPACAM2",
                    "GNG8","OLIG2","GNG13","RTP1","CXCR4", "SOX2","CYP2A13", "ERMN")
cell_type_names=c("Bowman’s glands",
                  "olfactory HBCs",
                  "olfactory ensheathing glia",
                  "olfactory microvillar cells",
                  "immature neurons",
                  "mature neurons",
                  "GBCs",
                  "sustentacular cells")

## Laoding paitent P2 and P3 matrices and formation of seurat object
batch_list=list("P2","P3")
batch_data_list=list("P2"=1,"P3"=1)
for( i in 1:length(batch_list))
{
  print(batch_list[[i]])
  s_object=Read10X(paste("~/corona_project/Input_files/",batch_list[[i]],sep=""))
  s_object=CreateSeuratObject(counts =s_object, min.cells = 0, min.features = 400, project = "P23")
  s_object[["percent.mt"]] <- PercentageFeatureSet(s_object, pattern = "^MT-")
  s_object <- subset(s_object, subset = nFeature_RNA >100 & nFeature_RNA <8000 & percent.mt <10)
  s_object@meta.data[, "run"] <- batch_list[i]
  s_object=NormalizeData(s_object)
  batch_data_list[[i]]=FindVariableFeatures(s_object, selection.method = "vst", nfeatures =5000)
}

## [1] "P2"
## [1] "P3"

batch_data_list

## $P2
## An object of class Seurat 
## 33538 features across 10881 samples within 1 assay 
## Active assay: RNA (33538 features, 5000 variable features)
## 
## $P3
## An object of class Seurat 
## 33538 features across 5415 samples within 1 assay 
## Active assay: RNA (33538 features, 5000 variable features)

saveRDS(batch_data_list,"~/corona_project/batch_data_list_5000_23.rds")

## Intigration and batch effect removing using Seurat
set.seed(1)
run.anchors <- FindIntegrationAnchors(object.list = batch_data_list, dims = 1:30,anchor.features = 5000)

## Computing 5000 integration features

## Scaling features for provided objects

## Finding all pairwise anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 11489 anchors

## Filtering anchors

##  Retained 8787 anchors

## Extracting within-dataset neighbors

a=run.anchors
run.combined <- IntegrateData(anchorset = run.anchors, dims = 1:30)

## Merging dataset 2 into 1

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Warning: Adding a command log without an assay associated with it

b=run.combined
sce_3_1_1_before=run.combined
saveRDS(sce_3_1_1_before,"~/corona_project/sce_3_1_1_before_5000_23.rds")

## Dim reduction and clustering
DefaultAssay(run.combined) <- "integrated"
run.combined <- ScaleData(run.combined, verbose = FALSE)
run.combined <- RunPCA(run.combined, npcs = 30, verbose = FALSE)
run.combined <- RunUMAP(run.combined, reduction = "pca", dims = 1:30)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 23:48:30 UMAP embedding parameters a = 0.9922 b = 1.112

## 23:48:30 Read 16296 rows and found 30 numeric columns

## 23:48:30 Using Annoy for neighbor search, n_neighbors = 30

## 23:48:30 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 23:48:33 Writing NN index file to temp file /tmp/RtmpQ1SbcN/file18d38f17eda
## 23:48:33 Searching Annoy index using 1 thread, search_k = 3000
## 23:48:37 Annoy recall = 100%
## 23:48:37 Commencing smooth kNN distance calibration using 1 thread
## 23:48:38 Initializing from normalized Laplacian + noise
## 23:48:39 Commencing optimization for 200 epochs, with 707164 positive edges
## 23:48:44 Optimization finished

run.combined <- FindNeighbors(run.combined, reduction = "pca", dims = 1:30)

## Computing nearest neighbor graph
## Computing SNN

run.combined <- FindClusters(run.combined, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 16296
## Number of edges: 654023
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9526
## Number of communities: 22
## Elapsed time: 2 seconds

## Cell-types annotation
for(i in 0:as.integer(length(marker_genes)/4))
{
  a=(i*4+1)
  b=((i+1)*4)
  if(i==as.integer(length(marker_genes)/4))
    b=length(marker_genes)
  if(a>length(marker_genes))
     next
  print(i)
  print(VlnPlot(run.combined, features = marker_genes[a:b],ncol=2))
  print(FeaturePlot(run.combined, features = marker_genes[a:b],cols=c("lightgrey", "red"),label = TRUE,pt.size=2,order = TRUE))
}

## [1] 0

## [1] 1

## [1] 2

## [1] 3

## [1] 4

## Warning: Could not find FOXL1 in the default search locations, found in RNA
## assay instead

## Warning: Could not find FOXL1 in the default search locations, found in RNA
## assay instead

## [1] 5

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## [1] 6

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## [1] 7

## [1] 8

## [1] 9

run.combined<-RenameIdents(run.combined,`13`="Immature neurons", `14`="Mature neurons", `3`="Olfactory HBCs",`19`="Olfactory microvillar cells",
                   `5`="Bowman's gland", `16`="Olfactory ensheathing glia", `9`="Sustentacular cells",`20`="GBCs" )
run.combined=subset(run.combined, idents = c("Immature neurons","Mature neurons","Olfactory HBCs","Olfactory microvillar cells","Bowman's gland",
                          "Olfactory ensheathing glia","Sustentacular cells","GBCs"))
sce_3_1_1_after=run.combined
saveRDS(sce_3_1_1_after,"~/corona_project/sce_3_1_1_after_5000_23.rds")

## Ploting heat map for selected gene markers on all clusters
sce_3_1_1_after_5000_23 <- readRDS("~/corona_project/sce_3_1_1_after_5000_23.rds")
cluster.averages <- AverageExpression(sce_3_1_1_after_5000_23, return.seurat = TRUE)

## Finished averaging RNA for cluster Immature neurons
## Finished averaging RNA for cluster Mature neurons
## Finished averaging RNA for cluster Olfactory HBCs
## Finished averaging RNA for cluster Olfactory microvillar cells
## Finished averaging RNA for cluster Bowman's gland
## Finished averaging RNA for cluster Olfactory ensheathing glia
## Finished averaging RNA for cluster Sustentacular cells
## Finished averaging RNA for cluster GBCs
## Finished averaging integrated for cluster Immature neurons
## Finished averaging integrated for cluster Mature neurons
## Finished averaging integrated for cluster Olfactory HBCs
## Finished averaging integrated for cluster Olfactory microvillar cells
## Finished averaging integrated for cluster Bowman's gland
## Finished averaging integrated for cluster Olfactory ensheathing glia
## Finished averaging integrated for cluster Sustentacular cells
## Finished averaging integrated for cluster GBCs
## Centering and scaling data matrix

DoHeatmap(sce_3_1_1_after_5000_23,features =c(marker_genes_temp),size = 3, draw.lines = FALSE )

DoHeatmap(cluster.averages,features =c(marker_genes_temp),size = 3, draw.lines = FALSE )

## Vln plots with all selected gene markers
plot=list()
for(i in 1:length(marker_genes_temp))
{
  print(i)
  if(i %in% c(1,5,9))
    plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
      NoLegend() +   theme(axis.title.x=element_blank(),
                           axis.text.x=element_blank(),
                           axis.ticks.x=element_blank(),
                           axis.line.x = element_blank())
  if(i %in% c(14,15,16))
    plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
      NoLegend() +   theme(axis.title.y=element_blank(),
                           axis.text.y=element_blank(),
                           axis.ticks.y=element_blank(),
                           axis.line.y = element_blank())
  if(i %in% c(13))
    plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
      NoLegend()
  if(i %in% c(2,3,4,6,7,8,10,11,12))
    plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
      NoLegend()+   theme(axis.title.x=element_blank(),
                          axis.text.x=element_blank(),
                          axis.ticks.x=element_blank(),
                          axis.line.x = element_blank(),
                          axis.title.y=element_blank(),
                          axis.text.y=element_blank(),
                          axis.ticks.y=element_blank(),
                          axis.line.y = element_blank())
  
}

## [1] 1

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [1] 6

## [1] 7

## [1] 8

## [1] 9

## [1] 10

## [1] 11

## [1] 12

## [1] 13

## [1] 14

## [1] 15

## [1] 16

cowplot::plot_grid(plotlist = plot,nrow=4,rel_heights = c(1,1,1,2))

## Feature plots with all selected gene markers
plot=list()
plot[[1]]=DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=TRUE)+
  labs(title = "") + NoAxes() + NoLegend()
for(i in 1:length(marker_genes_temp))
{
  print(i)
  plot[[i+1]]<-print(FeaturePlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]],cols=c("lightgrey", "red"),
                                 pt.size = point_size,label.size = label_size,order = TRUE)) +
    NoLegend() + NoAxes()
}

## [1] 1

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [1] 6

## [1] 7

## [1] 8

## [1] 9

## [1] 10

## [1] 11

## [1] 12

## [1] 13

## [1] 14

## [1] 15

## [1] 16

AT=cowplot::plot_grid(plotlist = plot[c(3,4,5)],ncol=3)
AB=cowplot::plot_grid(plotlist = plot[c(12,13,11)],ncol=3)
AR=cowplot::plot_grid(plotlist = plot[c(6,7,8,9,10)],nrow=5)
AL=cowplot::plot_grid(plotlist = plot[c(2,14,15,16,17)],nrow=5)
AC=cowplot::plot_grid(plotlist = plot[c(1)],nrow=1)
ATCB=cowplot::plot_grid(plotlist = list(AT,AC,AB),nrow=3,rel_heights = c(1,4,1))
cowplot::plot_grid(plotlist = list(AL,ATCB,AR),ncol=3,rel_widths = c(1,2,1))

### Dim plot to show all cell types with different colors
DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=TRUE)+ labs(title = "") + NoLegend()

DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=FALSE)+ labs(title = "")

## FeaturePlot to show expression of gene markers in respective cell markers
### FeaturePlot to show the expression of ACE2 in each cell type
FeaturePlot(run.combined, features = "ACE2",combine = FALSE,cols=c("lightgrey", "red"),pt.size = point_size,label.size = label_size,order =TRUE)

## Warning: Could not find ACE2 in the default search locations, found in RNA assay
## instead

## [[1]]

### FeaturePlot to show the expression of TMPRSS2 in each cell type
FeaturePlot(run.combined, features = "TMPRSS2",combine = FALSE,cols=c("lightgrey", "green"),pt.size = point_size,label.size = label_size,order =TRUE)

## [[1]]

### FeaturePlot to show the expression of CTSL in each cell type
FeaturePlot(run.combined, features = "CTSL",combine = FALSE,cols=c("lightgrey", "blue"),pt.size = point_size,label.size = label_size,order =TRUE)

## Warning: Could not find CTSL in the default search locations, found in RNA assay
## instead

## [[1]]

### FeaturePlot to show the expression of BSG in each cell type
FeaturePlot(run.combined, features = "BSG",combine = FALSE,cols=c("lightgrey", "magenta"),pt.size = point_size,label.size = label_size,order =TRUE)

## [[1]]

## Initiallizing a matrix for 24 * 8, 8 for cells type and 24 for different cases as shown below
count_matrix=matrix(0,24,8)
rownames(count_matrix)<-c("ACE2_count","ACE2_percent","TMPRSS2_count","TMPRSS2_percent",
                          "CTSL_count","CTSL_percent","BSG_count","BSG_percent",
                          "intersect_of_ACE2_TMPRSS2_count","ACE2_contri_percent","TMPRSS2_contri_percent",
                          "intersect_of_ACE2_CTSL_count","ACE2_contri_percent","CTSL_contri_percent",
                          "intersect_of_BSG_TMPRSS2_count","BSG_contri_percent","TMPRSS2_contri_percent",
                          "intersect_of_BSG_CTSL_count","BSG_contri_percent","CSTL_contri_percent",
                          "ACE2_TMPRSS2_CTSL","ACE2_TMPRSS2_CTSB",
                          "BSB_TMPRSS2_CTSL","BSG_TMPRSS2_CTSB")
colnames(count_matrix)<-unique(Idents(run.combined))
### Percentage of ACE2 expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[1,names(a)]=a
print(sum(a))

## [1] 32

a=a*100/sum(a)
count_matrix[2,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)

## [1] 33538  3906

### Percentage of TMPRSS2 expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
count_matrix[3,names(a)]=a
print(sum(a))

## [1] 666

a=a*100/sum(a)
count_matrix[4,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)

## [1] 33538  3906

### Percentage of CTSL expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
count_matrix[5,names(a)]=a
print(sum(a))

## [1] 1631

a=a*100/sum(a)
count_matrix[6,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)

## [1] 33538  3906

### Percentage of BSG expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[7,names(a)]=a
print(sum(a))

## [1] 3014

a=a*100/sum(a)
count_matrix[8,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)

## [1] 33538  3906

### Bar plot of intersect of both ACE2 and TMPRSS2
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[9,names(intersect_cells)]=intersect_cells
TMPRSS2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
ACE2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[10,names(ACE2_cells)]=(intersect_cells/ACE2_cells)*100
count_matrix[11,names(TMPRSS2_cells)]=(intersect_cells/TMPRSS2_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/TMPRSS2_cells)*100),((intersect_cells/ACE2_cells)*100)),
                    genes=c(rep("TMPRSS2",length(TMPRSS2_cells)),rep("ACE2",length(ACE2_cells))),
                    cell_types=c(names(TMPRSS2_cells),names(ACE2_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

## Warning: Removed 4 rows containing missing values (position_stack).

### Bar plot of intersect of both ACE2 and CTSL
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[12,names(intersect_cells)]=intersect_cells
CTSL_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
ACE2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[13,names(ACE2_cells)]=(intersect_cells/ACE2_cells)*100
count_matrix[14,names(CTSL_cells)]=(intersect_cells/CTSL_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/CTSL_cells)*100),((intersect_cells/ACE2_cells)*100)),
                    genes=c(rep("CTSL",length(CTSL_cells)),rep("ACE2",length(ACE2_cells))),
                    cell_types=c(names(CTSL_cells),names(ACE2_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

## Warning: Removed 4 rows containing missing values (position_stack).

### Bar plot of intersect of both BSG and TMPRSS2
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[15,names(intersect_cells)]=intersect_cells
TMPRSS2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
BSG_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[16,names(BSG_cells)]=(intersect_cells/BSG_cells)*100
count_matrix[17,names(TMPRSS2_cells)]=(intersect_cells/TMPRSS2_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/TMPRSS2_cells)*100),((intersect_cells/BSG_cells)*100)),
                    genes=c(rep("TMPRSS2",length(TMPRSS2_cells)),rep("BSG",length(BSG_cells))),
                    cell_types=c(names(TMPRSS2_cells),names(BSG_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar plot of intersect of both BSG and CTSL
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[18,names(intersect_cells)]=intersect_cells
CTSL_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
BSG_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[19,names(BSG_cells)]=(intersect_cells/BSG_cells)*100
count_matrix[20,names(CTSL_cells)]=(intersect_cells/CTSL_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/CTSL_cells)*100),((intersect_cells/BSG_cells)*100)),
                    genes=c(rep("CTSL",length(CTSL_cells)),rep("BSG",length(BSG_cells))),
                    cell_types=c(names(CTSL_cells),names(BSG_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
  geom_bar(stat="identity") +  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Intersect of ACE2, TMPRSS2 and CTSL
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[21,names(intersect_cells)]=intersect_cells
### Intersect of ACE2, TMPRSS2 and CTSB
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
col_name_CTSB=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSB",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSB))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[22,names(intersect_cells)]=intersect_cells
### Intersect of BSG, TMPRSS2 and CTSL
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[23,names(intersect_cells)]=intersect_cells
### Intersect of BSG, TMPRSS2 and CTSB
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
col_name_CTSB=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSB",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSB))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[24,names(intersect_cells)]=intersect_cells
count_matrix

##                                 Immature neurons Sustentacular cells
## ACE2_count                             0.0000000           15.000000
## ACE2_percent                           0.0000000           46.875000
## TMPRSS2_count                          5.0000000          345.000000
## TMPRSS2_percent                        0.7507508           51.801802
## CTSL_count                           167.0000000          360.000000
## CTSL_percent                          10.2391171           22.072348
## BSG_count                            347.0000000          622.000000
## BSG_percent                           11.5129396           20.637027
## intersect_of_ACE2_TMPRSS2_count        0.0000000            9.000000
## ACE2_contri_percent                          NaN           60.000000
## TMPRSS2_contri_percent                 0.0000000            2.608696
## intersect_of_ACE2_CTSL_count           0.0000000           11.000000
## ACE2_contri_percent                          NaN           73.333333
## CTSL_contri_percent                    0.0000000            3.055556
## intersect_of_BSG_TMPRSS2_count         5.0000000          334.000000
## BSG_contri_percent                     1.4409222           53.697749
## TMPRSS2_contri_percent               100.0000000           96.811594
## intersect_of_BSG_CTSL_count          148.0000000          347.000000
## BSG_contri_percent                    42.6512968           55.787781
## CSTL_contri_percent                   88.6227545           96.388889
## ACE2_TMPRSS2_CTSL                      0.0000000            7.000000
## ACE2_TMPRSS2_CTSB                      0.0000000            8.000000
## BSB_TMPRSS2_CTSL                       2.0000000          217.000000
## BSG_TMPRSS2_CTSB                       3.0000000          310.000000
##                                 Olfactory HBCs Mature neurons
## ACE2_count                          11.0000000      0.0000000
## ACE2_percent                        34.3750000      0.0000000
## TMPRSS2_count                       14.0000000      2.0000000
## TMPRSS2_percent                      2.1021021      0.3003003
## CTSL_count                         507.0000000    112.0000000
## CTSL_percent                        31.0852238      6.8669528
## BSG_count                          687.0000000    249.0000000
## BSG_percent                         22.7936297      8.2614466
## intersect_of_ACE2_TMPRSS2_count      0.0000000      0.0000000
## ACE2_contri_percent                  0.0000000            NaN
## TMPRSS2_contri_percent               0.0000000      0.0000000
## intersect_of_ACE2_CTSL_count         4.0000000      0.0000000
## ACE2_contri_percent                 36.3636364            NaN
## CTSL_contri_percent                  0.7889546      0.0000000
## intersect_of_BSG_TMPRSS2_count      13.0000000      2.0000000
## BSG_contri_percent                   1.8922853      0.8032129
## TMPRSS2_contri_percent              92.8571429    100.0000000
## intersect_of_BSG_CTSL_count        342.0000000    109.0000000
## BSG_contri_percent                  49.7816594     43.7751004
## CSTL_contri_percent                 67.4556213     97.3214286
## ACE2_TMPRSS2_CTSL                    0.0000000      0.0000000
## ACE2_TMPRSS2_CTSB                    0.0000000      0.0000000
## BSB_TMPRSS2_CTSL                     7.0000000      2.0000000
## BSG_TMPRSS2_CTSB                     6.0000000      1.0000000
##                                 Olfactory ensheathing glia Bowman's gland
## ACE2_count                                       0.0000000      4.0000000
## ACE2_percent                                     0.0000000     12.5000000
## TMPRSS2_count                                    6.0000000    239.0000000
## TMPRSS2_percent                                  0.9009009     35.8858859
## CTSL_count                                      90.0000000    281.0000000
## CTSL_percent                                     5.5180871     17.2286941
## BSG_count                                      157.0000000    772.0000000
## BSG_percent                                      5.2090246     25.6138023
## intersect_of_ACE2_TMPRSS2_count                  0.0000000      1.0000000
## ACE2_contri_percent                                    NaN     25.0000000
## TMPRSS2_contri_percent                           0.0000000      0.4184100
## intersect_of_ACE2_CTSL_count                     0.0000000      2.0000000
## ACE2_contri_percent                                    NaN     50.0000000
## CTSL_contri_percent                              0.0000000      0.7117438
## intersect_of_BSG_TMPRSS2_count                   5.0000000    228.0000000
## BSG_contri_percent                               3.1847134     29.5336788
## TMPRSS2_contri_percent                          83.3333333     95.3974895
## intersect_of_BSG_CTSL_count                     68.0000000    266.0000000
## BSG_contri_percent                              43.3121019     34.4559585
## CSTL_contri_percent                             75.5555556     94.6619217
## ACE2_TMPRSS2_CTSL                                0.0000000      0.0000000
## ACE2_TMPRSS2_CTSB                                0.0000000      1.0000000
## BSB_TMPRSS2_CTSL                                 2.0000000     91.0000000
## BSG_TMPRSS2_CTSB                                 1.0000000    215.0000000
##                                 Olfactory microvillar cells       GBCs
## ACE2_count                                         0.000000   2.000000
## ACE2_percent                                       0.000000   6.250000
## TMPRSS2_count                                     16.000000  39.000000
## TMPRSS2_percent                                    2.402402   5.855856
## CTSL_count                                        59.000000  55.000000
## CTSL_percent                                       3.617413   3.372164
## BSG_count                                        113.000000  67.000000
## BSG_percent                                        3.749171   2.222960
## intersect_of_ACE2_TMPRSS2_count                    0.000000   2.000000
## ACE2_contri_percent                                     NaN 100.000000
## TMPRSS2_contri_percent                             0.000000   5.128205
## intersect_of_ACE2_CTSL_count                       0.000000   1.000000
## ACE2_contri_percent                                     NaN  50.000000
## CTSL_contri_percent                                0.000000   1.818182
## intersect_of_BSG_TMPRSS2_count                    16.000000  37.000000
## BSG_contri_percent                                14.159292  55.223881
## TMPRSS2_contri_percent                           100.000000  94.871795
## intersect_of_BSG_CTSL_count                       57.000000  50.000000
## BSG_contri_percent                                50.442478  74.626866
## CSTL_contri_percent                               96.610169  90.909091
## ACE2_TMPRSS2_CTSL                                  0.000000   1.000000
## ACE2_TMPRSS2_CTSB                                  0.000000   1.000000
## BSB_TMPRSS2_CTSL                                   9.000000  32.000000
## BSG_TMPRSS2_CTSB                                  11.000000  33.000000

## Bar plots to show + expressed cells of some gene markers as shown beelow
### Bar side plot of ACE2, TMPRSS2 and (ACE2 and TMPRSS2)
count=c(count_matrix["ACE2_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_ACE2_TMPRSS2_count",])
count[count==0]<-1
A_T_AT_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("ACE2",8),rep("TMPRSS2",8),rep("Intersect (ACE2 and TMPRSS2)",8)))
ggplot(data=A_T_AT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity", position=position_dodge())+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["BSG_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_BSG_CTSL_count",])
count[count==0]<-1
B_C_BC_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("BSG",8),rep("CTSL",8),rep("Intersect (BSG and CTSL)",8)))
ggplot(data=B_C_BC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity", position=position_dodge())+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["ACE2_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_ACE2_CTSL_count",])
count[count==0]<-1
A_C_AC_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("ACE2",8),rep("CTSL",8),rep("Intersect (ACE2 and CTSL)",8)))
ggplot(data=A_C_AC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity", position=position_dodge())+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, TMPRSS2 and (BSG and TMPRSS2)
count=c(count_matrix["BSG_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_BSG_TMPRSS2_count",])
count[count==0]<-1
B_T_BT_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("BSG",8),rep("TMPRSS2",8),rep("Intersect (BSG and TMPRSS2)",8)))
ggplot(data=B_T_BT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity", position=position_dodge())+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of ACE2, TMPRSS2 and (ACE2 and TMPRSS2)
count=c(count_matrix["ACE2_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_ACE2_TMPRSS2_count",])
count[count==0]<-1
A_T_AT_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("ACE2",8),rep("TMPRSS2",8),rep("Intersect (ACE2 and TMPRSS2)",8)))
ggplot(data=A_T_AT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity")+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["BSG_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_BSG_CTSL_count",])
count[count==0]<-1
B_C_BC_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("BSG",8),rep("CTSL",8),rep("Intersect (BSG and CTSL)",8)))
ggplot(data=B_C_BC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity")+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["ACE2_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_ACE2_CTSL_count",])
count[count==0]<-1
A_C_AC_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("ACE2",8),rep("CTSL",8),rep("Intersect (ACE2 and CTSL)",8)))
ggplot(data=A_C_AC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity")+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, TMPRSS2 and (BSG and TMPRSS2)
count=c(count_matrix["BSG_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_BSG_TMPRSS2_count",])
count[count==0]<-1
B_T_BT_df=data.frame("count"=count,"cell_markers"= names(count),
                     "gene_markers"=c(rep("BSG",8),rep("TMPRSS2",8),rep("Intersect (BSG and TMPRSS2)",8)))
ggplot(data=B_T_BT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
  geom_bar(stat="identity")+
  theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
  scale_y_continuous(trans='log10') +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

## DE genes analysis for selected cell types
run.combined=sce_3_1_1_after_5000_23
DefaultAssay(run.combined) <- "RNA"
seurat_subset_list=list()
p_val=0.05
l_fold=1
### DE genes analysis on ACE2+ and ACE2-
gene_marker="ACE2"
print(sort(count_matrix[1,]))

##            Immature neurons              Mature neurons 
##                           0                           0 
##  Olfactory ensheathing glia Olfactory microvillar cells 
##                           0                           0 
##                        GBCs              Bowman's gland 
##                           2                           4 
##              Olfactory HBCs         Sustentacular cells 
##                          11                          15

selected_cells_type=sort(count_matrix[1,])
selected_cells_type=selected_cells_type[selected_cells_type>2]
de_methods=c("wilcox","MAST","bimod","t","negbinom","poisson","LR")
for(k in de_methods)
  for(i in 1:length(selected_cells_type))
  {
    print(paste(k,"_",i,sep=""))
    seurat_subset=subset(run.combined, idents = names(selected_cells_type)[i])
    positive_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker,])>0)])
    negative_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker,])<=0)])
    pos_cells=paste(gene_marker,"+",sep="")
    neg_cells=paste(gene_marker,"-",sep="")
    pos=c(rep(pos_cells,length(positive_cells_name)))
    names(pos)=positive_cells_name
    neg=c(rep(neg_cells,length(negative_cells_name)))
    names(neg)=negative_cells_name
    pos_neg=c(pos,neg)
    print(paste("cells in seurat object after subsetting with ",names(selected_cells_type)[i]," = ",dim(seurat_subset)[2],sep=""))
    print(paste(names(selected_cells_type)[i]," total cells = ",sum(pos_neg[which(names(pos_neg)%in%rownames(seurat_subset@meta.data))]==pos_neg),sep=""))
    print(paste(pos_cells," cells = ",sum(pos_neg==pos_cells),sep=""))
    print(paste(neg_cells," cells = ",sum(pos_neg==neg_cells),sep=""))
    seurat_subset=AddMetaData(
      object = seurat_subset,
      metadata = pos_neg,
      col.name = 'seurat_subset'
    )
    seurat_subset_list[[paste(names(selected_cells_type)[i],"_",gene_marker,sep="")]]<-seurat_subset
    Idents(seurat_subset)<-pos_neg
    seurat_subset_genes=FindMarkers(seurat_subset, ident.1=pos_cells,ident.2=neg_cells,test.use = k,min.cells.group=1)
    print(seurat_subset_genes[1:2,])
    seurat_subset_genes_up=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC>l_fold & seurat_subset_genes$p_val<p_val]
    seurat_subset_genes_down=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC<(-l_fold) & seurat_subset_genes$p_val<p_val]
    seurat_subset_genes_names=c(seurat_subset_genes_up,seurat_subset_genes_down)
    print(paste(names(selected_cells_type)[i],"_",pos_cells," genes = ",length(seurat_subset_genes_up),
                " and ",names(selected_cells_type)[i],"_",neg_cells," genes = ",length(seurat_subset_genes_down),sep=""))
    DE_info=data.frame("log2FC"=as.numeric(seurat_subset_genes$avg_logFC),"p_val"=as.numeric(seurat_subset_genes$p_val))
    lab=rownames(seurat_subset_genes)
    print(EnhancedVolcano(DE_info,
                          lab = lab,
                          x = 'log2FC',
                          y = 'p_val',
                          xlim = c(-2, 2),
                          title = names(selected_cells_type)[i],
                          pCutoff = p_val,
                          FCcutoff = l_fold,
                          colAlpha = 1))
  }

## [1] "wilcox_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
##                   p_val avg_logFC pct.1 pct.2    p_val_adj
## AP003086.2 1.069865e-50 0.3138844  0.25 0.000 3.588112e-46
## AC009163.1 4.596423e-26 0.3134687  0.25 0.001 1.541548e-21
## [1] "Bowman's gland_ACE2+ genes = 3 and Bowman's gland_ACE2- genes = 2"

## [1] "wilcox_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
##               p_val avg_logFC pct.1 pct.2    p_val_adj
## SCN11A 1.905727e-15 0.3122231 0.182 0.004 6.391426e-11
## CHODL  1.521107e-07 0.2897168 0.182 0.010 5.101488e-03
## [1] "Olfactory HBCs_ACE2+ genes = 0 and Olfactory HBCs_ACE2- genes = 2"

## [1] "wilcox_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
##               p_val avg_logFC pct.1 pct.2 p_val_adj
## CARTPT 1.257096e-05 0.2589717   0.2 0.022 0.4216048
## MYCL   2.385754e-05 0.6544713   0.8 0.303 0.8001343
## [1] "Sustentacular cells_ACE2+ genes = 0 and Sustentacular cells_ACE2- genes = 0"

## [1] "MAST_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"

## Assuming data assay in position 1, with name et is log-transformed.
## 
## Done!
## Combining coefficients and standard errors

## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.

## Calculating log-fold changes

## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.

## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!

## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

##                  p_val avg_logFC pct.1 pct.2 p_val_adj
## R3HDM4    0.0007098688 0.7573622   1.0 0.205         1
## LINC01752 0.0008895931 0.4610704   0.5 0.007         1
## [1] "Bowman's gland_ACE2+ genes = 5 and Bowman's gland_ACE2- genes = 2"

## [1] "MAST_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"

## Assuming data assay in position 1, with name et is log-transformed.
## 
## Done!
## Combining coefficients and standard errors

## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.

## Calculating log-fold changes

## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.

## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!

## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

##                p_val avg_logFC pct.1 pct.2 p_val_adj
## CTSA    9.283015e-06 0.4894748 0.182 0.243 0.3113337
## SMARCC2 3.723033e-05 0.5442944 0.273 0.235 1.0000000
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 2"

## [1] "MAST_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"

## Assuming data assay in position 1, with name et is log-transformed.
## 
## Done!
## Combining coefficients and standard errors

## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.

## Calculating log-fold changes

## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.

## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!

## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

##                   p_val avg_logFC pct.1 pct.2    p_val_adj
## AC009113.1 3.740250e-11 0.5956145 0.067 0.104 1.254405e-06
## C8orf58    4.486527e-05 0.5854481 0.200 0.112 1.000000e+00
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 0"

## [1] "bimod_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
##                 p_val avg_logFC pct.1 pct.2   p_val_adj
## C19orf70 1.182295e-07 0.5201403  0.75 0.621 0.003965182
## PPP1R32  6.583710e-07 0.3598497  0.50 0.012 0.022080447
## [1] "Bowman's gland_ACE2+ genes = 2 and Bowman's gland_ACE2- genes = 3"

## [1] "bimod_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
##              p_val avg_logFC pct.1 pct.2    p_val_adj
## IGLC3 6.825386e-13  2.190882 0.364 0.205 2.289098e-08
## IGLC2 5.545072e-10  1.824826 0.364 0.398 1.859706e-05
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 2"

## [1] "bimod_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
##                   p_val avg_logFC pct.1 pct.2    p_val_adj
## URB1-AS1   2.788252e-16 0.8820008 0.600 0.269 9.351240e-12
## AL034549.1 2.358246e-15 0.6249744 0.267 0.118 7.909085e-11
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 2"

## [1] "t_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
##                  p_val  avg_logFC pct.1 pct.2     p_val_adj
## HDLBP    1.724021e-155 -0.8721947     0 0.638 5.782021e-151
## MPHOSPH8 1.049661e-153 -0.8076480     0 0.628 3.520352e-149
## [1] "Bowman's gland_ACE2+ genes = 0 and Bowman's gland_ACE2- genes = 9"

## [1] "t_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
##                 p_val  avg_logFC pct.1 pct.2     p_val_adj
## LAMB1   9.368870e-135 -0.7778893     0 0.486 3.142132e-130
## ATP13A3 1.615811e-103 -0.5378039     0 0.393  5.419108e-99
## [1] "Olfactory HBCs_ACE2+ genes = 0 and Olfactory HBCs_ACE2- genes = 3"

## [1] "t_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
##                p_val  avg_logFC pct.1 pct.2    p_val_adj
## TMEM238 8.163443e-25 -0.3774210     0 0.190 2.737855e-20
## XIST    1.701171e-17 -0.5124281     0 0.109 5.705389e-13
## [1] "Sustentacular cells_ACE2+ genes = 0 and Sustentacular cells_ACE2- genes = 2"

## [1] "negbinom_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in FUN(X[[i]], ...): Skipping gene --- 649. Fewer than 3 cells in both
## clusters.

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

##                  p_val avg_logFC pct.1 pct.2   p_val_adj
## LINC01752 1.263145e-07 0.3988132   0.5 0.007 0.004236336
## PPP1R32   1.402014e-06 0.3933035   0.5 0.012 0.047020752
## [1] "Bowman's gland_ACE2+ genes = 0 and Bowman's gland_ACE2- genes = 7"

## [1] "negbinom_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached

##               p_val avg_logFC pct.1 pct.2    p_val_adj
## IGLC3  1.960891e-10  1.797373 0.364 0.205 6.576436e-06
## MT-ND1 9.184367e-05 -1.105728 1.000 0.995 1.000000e+00
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 9"

## [1] "negbinom_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

##               p_val avg_logFC pct.1 pct.2   p_val_adj
## BICDL1 6.437705e-08 0.3978762 0.467 0.071 0.002159077
## MRPL48 9.001799e-07 0.6287818 0.733 0.355 0.030190234
## [1] "Sustentacular cells_ACE2+ genes = 5 and Sustentacular cells_ACE2- genes = 2"

## [1] "poisson_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"

## Warning in FUN(X[[i]], ...): Skipping gene --- 649. Fewer than 3 cells in both
## clusters.

##              p_val avg_logFC pct.1 pct.2     p_val_adj
## SAA1 7.245093e-219  2.152606  0.25 0.271 2.429859e-214
## LYZ  3.583759e-134 -3.023938  0.75 0.934 1.201921e-129
## [1] "Bowman's gland_ACE2+ genes = 14 and Bowman's gland_ACE2- genes = 23"

## [1] "poisson_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
##               p_val avg_logFC pct.1 pct.2     p_val_adj
## IGLC3 1.832406e-132  1.797373 0.364 0.205 6.145523e-128
## IGLC2  1.357811e-65  1.425866 0.364 0.398  4.553827e-61
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 9"

## [1] "poisson_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
##                p_val avg_logFC pct.1 pct.2     p_val_adj
## BPIFB1  0.000000e+00 2.5050723 0.133 0.062  0.000000e+00
## GSTP1  4.435057e-238 0.4987906 1.000 0.992 1.487429e-233
## [1] "Sustentacular cells_ACE2+ genes = 14 and Sustentacular cells_ACE2- genes = 2"

## Warning: One or more p-values is 0. Converting to 10^-1 * current lowest non-
## zero p-value...

## [1] "LR_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                  p_val avg_logFC pct.1 pct.2 p_val_adj
## LINC01752 6.320178e-05 0.4610704  0.50 0.007         1
## C3orf33   3.938426e-04 0.5344172  0.75 0.053         1
## [1] "Bowman's gland_ACE2+ genes = 4 and Bowman's gland_ACE2- genes = 3"

## [1] "LR_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##              p_val  avg_logFC pct.1 pct.2 p_val_adj
## TULP3 8.050212e-05  0.7332917 0.455 0.096         1
## ARL4D 1.345051e-04 -1.1427918 0.091 0.589         1
## [1] "Olfactory HBCs_ACE2+ genes = 1 and Olfactory HBCs_ACE2- genes = 3"

## [1] "LR_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                 p_val avg_logFC pct.1 pct.2 p_val_adj
## MYCL     0.0004320384 0.6544713   0.8 0.303         1
## URB1-AS1 0.0004526378 0.8820008   0.6 0.269         1
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 1"

### DE genes analysis on intersection of (ACE2 and TMPRSS2)+ and (ACE2 and TMPRSS2)-
gene_marker=c("ACE2","TMPRSS2")
print(sort(count_matrix[1,]))

##            Immature neurons              Mature neurons 
##                           0                           0 
##  Olfactory ensheathing glia Olfactory microvillar cells 
##                           0                           0 
##                        GBCs              Bowman's gland 
##                           2                           4 
##              Olfactory HBCs         Sustentacular cells 
##                          11                          15

selected_cells_type=sort(count_matrix[9,])
selected_cells_type=selected_cells_type[selected_cells_type>2]
for(k in de_methods)
  for(i in 1:length(selected_cells_type))
  {
    print(paste(k,"_",i,sep=""))
    seurat_subset=subset(run.combined, idents = names(selected_cells_type)[i])
    positive_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker[1],])>0 & as.matrix(seurat_subset@assays$RNA[gene_marker[2],])>0)])
    negative_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker[1],])<=0 | as.matrix(seurat_subset@assays$RNA[gene_marker[2],])<=0)])
    pos_cells=paste(gene_marker[1],"_and_",gene_marker[2],"+",sep="")
    neg_cells=paste(gene_marker[1],"_and_",gene_marker[2],"-",sep="")
    pos=c(rep(pos_cells,length(positive_cells_name)))
    names(pos)=positive_cells_name
    neg=c(rep(neg_cells,length(negative_cells_name)))
    names(neg)=negative_cells_name
    pos_neg=c(pos,neg)
    print(paste("cells in seurat object after subsetting with ",names(selected_cells_type)[i]," = ",dim(seurat_subset)[2],sep=""))
    print(paste(names(selected_cells_type)[i]," total cells = ",sum(pos_neg[which(names(pos_neg)%in%rownames(seurat_subset@meta.data))]==pos_neg),sep=""))
    print(paste(pos_cells," cells = ",sum(pos_neg==pos_cells),sep=""))
    print(paste(neg_cells," cells = ",sum(pos_neg==neg_cells),sep=""))
    seurat_subset=AddMetaData(
      object = seurat_subset,
      metadata = pos_neg,
      col.name = 'seurat_subset'
    )
    seurat_subset_list[[paste(names(selected_cells_type)[i],"_",gene_marker[1],"_and_",gene_marker[2],sep="")]]<-seurat_subset
    Idents(seurat_subset)<-pos_neg
    seurat_subset_genes=FindMarkers(seurat_subset, ident.1=pos_cells,ident.2=neg_cells,test.use = k,min.cells.group=1)
    print(seurat_subset_genes[1:2,])
    seurat_subset_genes_up=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC>l_fold & seurat_subset_genes$p_val<p_val]
    seurat_subset_genes_down=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC<(-l_fold) & seurat_subset_genes$p_val<p_val]
    seurat_subset_genes_names=c(seurat_subset_genes_up,seurat_subset_genes_down)
    print(paste(names(selected_cells_type)[i],"_",pos_cells," genes = ",length(seurat_subset_genes_up),
                " and ",names(selected_cells_type)[i],"_",neg_cells," genes = ",length(seurat_subset_genes_down),sep=""))
    DE_info=data.frame("log2FC"=as.numeric(seurat_subset_genes$avg_logFC),"p_val"=as.numeric(seurat_subset_genes$p_val))
    lab=rownames(seurat_subset_genes)
    print(EnhancedVolcano(DE_info,
                          lab = lab,
                          x = 'log2FC',
                          y = 'p_val',
                          xlim = c(-5, 5),
                          title = names(selected_cells_type)[i],
                          pCutoff = p_val,
                          FCcutoff = l_fold,
                          colAlpha = 1))
  }

## [1] "wilcox_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
##               p_val avg_logFC pct.1 pct.2  p_val_adj
## PLXNA1 1.104891e-06 0.2963436 0.778 0.169 0.03705582
## MTX3   2.878126e-06 0.3344437 0.778 0.201 0.09652658
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 0 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## [1] "MAST_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"

## Assuming data assay in position 1, with name et is log-transformed.
## 
## Done!
## Combining coefficients and standard errors

## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.

## Calculating log-fold changes

## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.

## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!

## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

##                 p_val avg_logFC pct.1 pct.2    p_val_adj
## C1orf159 1.818218e-09 0.7760622 0.111 0.081 6.097939e-05
## C8orf58  6.750627e-09 0.8492568 0.222 0.113 2.264025e-04
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 3 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## [1] "bimod_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
##                p_val avg_logFC pct.1 pct.2    p_val_adj
## SELENOO 5.270310e-15 0.7256710 0.333 0.209 1.767557e-10
## FAM19A4 2.976249e-12 0.7768982 0.222 0.098 9.981744e-08
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 5 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"

## [1] "t_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
##                p_val  avg_logFC pct.1 pct.2    p_val_adj
## MYL9    9.828220e-26 -0.2544036     0 0.227 3.296189e-21
## TMEM238 8.437425e-25 -0.3748267     0 0.188 2.829743e-20
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 0 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"

## [1] "negbinom_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in sqrt(1/i): NaNs produced

## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in FUN(X[[i]], ...): Skipping gene --- 1650. Fewer than 3 cells in both
## clusters.

## Warning in FUN(X[[i]], ...): Skipping gene --- 1780. Fewer than 3 cells in both
## clusters.

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in FUN(X[[i]], ...): Skipping gene --- 1957. Fewer than 3 cells in both
## clusters.

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

##               p_val avg_logFC pct.1 pct.2  p_val_adj
## KLHL32 4.679348e-07 0.2767530 0.333 0.011 0.01569360
## NBPF15 1.775853e-06 0.4278287 0.444 0.084 0.05955854
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 9 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## [1] "poisson_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"

## Warning in FUN(X[[i]], ...): Skipping gene --- 1650. Fewer than 3 cells in both
## clusters.

## Warning in FUN(X[[i]], ...): Skipping gene --- 1780. Fewer than 3 cells in both
## clusters.

## Warning in FUN(X[[i]], ...): Skipping gene --- 1957. Fewer than 3 cells in both
## clusters.

##         p_val avg_logFC pct.1 pct.2 p_val_adj
## SCGB1A1     0  3.269476 0.111 0.058         0
## BPIFB1      0  2.985891 0.111 0.063         0
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 27 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"

## Warning: One or more p-values is 0. Converting to 10^-1 * current lowest non-
## zero p-value...

## [1] "LR_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                p_val avg_logFC pct.1 pct.2 p_val_adj
## COPS7A  0.0004655404 0.7991358 0.889 0.438         1
## CCDC184 0.0004961425 0.3207393 0.333 0.052         1
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 3 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## Finding Stouffer score after log transfer and zscore
### Loading seaurat object
sce_3_1_1_after_5000_23 <- readRDS("~/corona_project/sce_3_1_1_after_5000_23.rds")
mat<-sce_3_1_1_after_5000_23@assays$RNA@counts
### Cell filtering
c_s <- Matrix::colSums(mat>=3)
l1 = stats::quantile(c_s,probs = 0.001)
l2 = stats::quantile(c_s,probs = 1)
keep_cells = intersect(which(c_s>=l1), which(c_s<=l2))
mat<-mat[,keep_cells]
### Gene filtering
r_s <- Matrix::rowSums(mat > 2)
keep_genes = which(r_s > 3)
mat<-mat[keep_genes,]
### Median normalization
r_s<-Matrix::rowSums(t(mat))
r_s_med<-stats::median(r_s)
mat<-Matrix::t(t(mat)/(r_s/r_s_med))
### Loading media genes file and intersecting genes and filtering matrix common genes
media_genes <- read_excel("~/corona_project/media-6.xlsx")$...3[-1]

## New names:
## * `` -> ...2
## * `` -> ...3
## * `` -> ...5
## * `` -> ...6
## * `` -> ...7
## * ...

common_genes=intersect(rownames(mat),media_genes)
mat=mat[common_genes,]
### Further insuring have cells with atleast 10% expressed genes
mat<-mat[,Matrix::colSums(mat>0)>(dim(mat)[1]/10)]
print(dim(mat))

## [1]  293 3493

seurat_RNA_mat<-mat
print(sum(apply(seurat_RNA_mat,1,function(x) sum(x)==0)))

## [1] 0

print(sum(apply(seurat_RNA_mat,2,function(x) sum(x)==0)))

## [1] 0

print(dim(seurat_RNA_mat))

## [1]  293 3493

### First psedo count 1 then log2 then zscore using scale function in R
seurat_RNA_mat<-scale(log2(seurat_RNA_mat+1),center = TRUE, scale = TRUE)
### Stouffer score using corto package
stouffer_score<- apply(seurat_RNA_mat,2,corto::stouffer)
print(sum(is.na(stouffer_score)))

## [1] 0

cell_types<-Idents(sce_3_1_1_after_5000_23)[names(stouffer_score)]
unique_cell_types=unique(cell_types)
Stouffer_score_df<-data.frame("Stouffer_score"=stouffer_score, "cell_types"=cell_types)
med=c()
med_n=c()
ct=unique(cell_types)
for(i in ct)
{
  med = c(med,median(stouffer_score[cell_types==i]))
  med_n = c(med_n,i)
  
}
med_df=data.frame("cell_types"=med_n, "med"=med)
ggplot(med_df, aes(x=cell_types, y=med, fill=cell_types)) +
  geom_bar(stat="identity")+theme_minimal() +
  theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))